rm(list=ls())
library(foreign)
library(dplyr)
library(ggplot2)
library(spatialEco)
library(scales)
library(gtable)
library(grid)
library(readstata13)

mydata <- read.dta("Drought_PC.dta")

Voteshare_inc_pm1 <- subset(mydata,  demean_prev_monsoon2 > -1 & demean_prev_monsoon2 < 1 & demean_inc_pm!="" & pm_align == 1)
Voteshare_inc_pm0 <- subset(mydata,  demean_prev_monsoon2 > -1 & demean_prev_monsoon2 < 1 & demean_inc_pm!="" & pm_align == 0)
Voteshare_inc_nonpm1 <- subset(mydata,  demean_prev_monsoon2 > -1 & demean_prev_monsoon2 < 1 & demean_inc_nonpm!="" & pm_align == 1)
Voteshare_inc_nonpm0 <- subset(mydata,  demean_prev_monsoon2 > -1 & demean_prev_monsoon2 < 1 & demean_inc_nonpm!="" & pm_align == 0)


############################
#Appendix Section A13 PM-MP Alignment
############################

############################
#Figure A18 (a)
############################

fit <- lm(demean_inc_pm ~ demean_prev_monsoon2 + I(demean_prev_monsoon2^2), data = Voteshare_inc_pm1)
prd <- data.frame(demean_prev_monsoon2 = seq(from = range(Voteshare_inc_pm1$demean_prev_monsoon2)[1], to = range(Voteshare_inc_pm1$demean_prev_monsoon2)[2], length.out = 100))
err <- predict(fit, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit

quad_Voteshare_inc_pm1<-ggplot(prd, aes(x = demean_prev_monsoon2, y = fit)) +
  theme_bw() +
  geom_line() +
  geom_smooth(aes(ymin = lci, ymax = uci), stat = "identity",colour = "black" ) +
  coord_fixed(ratio=0.1) +
  scale_y_continuous(limits = c(-9, 10),breaks = c(-8,-6,-4,-2,0,2,4,6,8,10) ) + 
  #geom_line(aes(y = lo.Voteshare_inc_pm$up.lim)) + geom_line(aes(y = lo.Voteshare_inc_pm$low.lim)) +
  geom_ribbon(data=prd, aes(ymin=lci ,ymax=uci), fill="grey50", alpha="0.5") +
  labs(y = "Vote Share (PM-Incumbent, demeaned)", x="Prev. Monsoon (demeaned)", title="Quadratic Prediction Plot: Vote Share (PM-Incumbent, PM-MP)") + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.title.y=element_text(vjust=-1),
        axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #panel.border = element_blank(),
        #panel.background = element_blank(),
        plot.margin = unit(c(0.3,0.7,0,0), "cm")) 

gg_Voteshare_inc_pm1<-ggplot(data=Voteshare_inc_pm1, aes(demean_prev_monsoon2)) +  
  geom_histogram(aes(y=..count../sum(..count..)), colour = "grey50", fill = "grey50", breaks=seq(-1, 1, by=0.2), alpha = "0.2") + 
  labs(y = "Share of Observations") +  
  coord_fixed(ratio=1) +
  scale_y_continuous(limits = c(0, 0.4), breaks = c(0,0.1,0.2,0.3,0.4)) + 
  theme(axis.title.y=element_text(vjust=-0.5),
        axis.text.y = element_text(vjust=0),
        #panel.background = element_rect(fill = NA),
        axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.ticks.length=unit(0.2,"cm"),
        plot.margin = unit(c(0.3,0.7,0,0), "cm")) 

grid.draw(quad_Voteshare_inc_pm1)

g1 <- ggplot_gtable(ggplot_build(quad_Voteshare_inc_pm1))
gpm1 <- ggplot_gtable(ggplot_build(quad_Voteshare_inc_pm1))
g2 <- ggplot_gtable(ggplot_build(gg_Voteshare_inc_pm1))

pp <- c(subset(g1$layout, name == "panel", se = t:r))
g <- gtable_add_grob(g1, g2$grobs[[which(g2$layout$name == "panel")]], pp$t, 
                     pp$l, pp$b, pp$l)

ia <- which(g2$layout$name == "axis-l")
ga <- g2$grobs[[ia]]
ax <- ga$children[[2]]
ax$widths <- rev(ax$widths)
ax$grobs <- rev(ax$grobs)
ax$grobs[[1]]$x <- ax$grobs[[1]]$x - unit(1, "npc") + unit(0.15, "cm")
g <- gtable_add_cols(g, g2$widths[g2$layout[ia, ]$l], length(g$widths) - 1)
g <- gtable_add_grob(g, ax, pp$t, length(g$widths) - 1, pp$b)

ia2 <- which(g2$layout$name == "ylab-l")
ga2 <- g2$grobs[[ia2]]
ga2$rot <- 90
g <- gtable_add_cols(g, g2$widths[g2$layout[ia2, ]$l], length(g$widths) - 1)
g <- gtable_add_grob(g, ga2, pp$t, length(g$widths) - 1, pp$b)

# draw it

grid.draw(g)

pdf("Voteshare_inc_pm1_total_boot_quad1_firststage_ms2.pdf") # Open a new pdf file
grid.draw(g) # Write the grid.arrange in the file
dev.off() # Close the file

############################
#Figure A18 (b)
############################


fit <- lm(demean_inc_pm ~ demean_prev_monsoon2 + I(demean_prev_monsoon2^2), data = Voteshare_inc_pm0)
prd <- data.frame(demean_prev_monsoon2 = seq(from = range(Voteshare_inc_pm0$demean_prev_monsoon2)[1], to = range(Voteshare_inc_pm0$demean_prev_monsoon2)[2], length.out = 100))
err <- predict(fit, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit

quad_Voteshare_inc_pm2<-ggplot(prd, aes(x = demean_prev_monsoon2, y = fit)) +
  theme_bw() +
  geom_line() +
  geom_smooth(aes(ymin = lci, ymax = uci), stat = "identity",colour = "black" ) +
  coord_fixed(ratio=0.1) +
  scale_y_continuous(limits = c(-9, 10),breaks = c(-8,-6,-4,-2,0,2,4,6,8,10) ) + 
  #geom_line(aes(y = lo.Voteshare_inc_pm$up.lim)) + geom_line(aes(y = lo.Voteshare_inc_pm$low.lim)) +
  geom_ribbon(data=prd, aes(ymin=lci ,ymax=uci), fill="grey50", alpha="0.5") +
  labs(y = "Vote Share (PM-Incumbent, demeaned)", x="Prev. Monsoon (demeaned)", title="Quadratic Prediction Plot: Vote Share (PM-Incumbent, non-PM-MP)") + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.title.y=element_text(vjust=-1),
        axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #panel.border = element_blank(),
        #panel.background = element_blank(),
        plot.margin = unit(c(0.3,0.7,0,0), "cm")) 

gg_Voteshare_inc_pm2<-ggplot(data=Voteshare_inc_pm0, aes(demean_prev_monsoon2)) +  
  geom_histogram(aes(y=..count../sum(..count..)), colour = "grey50", fill = "grey50", breaks=seq(-1, 1, by=0.2), alpha = "0.2") + 
  labs(y = "Share of Observations") +  
  coord_fixed(ratio=0.1) +
  scale_y_continuous(limits = c(0, 0.4), breaks = c(0,0.1,0.2,0.3,0.4)) + 
  theme(axis.title.y=element_text(vjust=-0.5),
        axis.text.y = element_text(vjust=0),
        panel.background = element_rect(fill = NA),
        axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        #panel.background = element_blank(),
        axis.ticks.length=unit(0.2,"cm"),
        plot.margin = unit(c(0.3,0.7,0,0), "cm")) 


g1 <- ggplot_gtable(ggplot_build(quad_Voteshare_inc_pm2))
gpm2 <- ggplot_gtable(ggplot_build(quad_Voteshare_inc_pm2))
g2 <- ggplot_gtable(ggplot_build(gg_Voteshare_inc_pm2))

pp <- c(subset(g1$layout, name == "panel", se = t:r))
g <- gtable_add_grob(g1, g2$grobs[[which(g2$layout$name == "panel")]], pp$t, 
                     pp$l, pp$b, pp$l)

ia <- which(g2$layout$name == "axis-l")
ga <- g2$grobs[[ia]]
ax <- ga$children[[2]]
ax$widths <- rev(ax$widths)
ax$grobs <- rev(ax$grobs)
ax$grobs[[1]]$x <- ax$grobs[[1]]$x - unit(1, "npc") + unit(0.15, "cm")
g <- gtable_add_cols(g, g2$widths[g2$layout[ia, ]$l], length(g$widths) - 1)
g <- gtable_add_grob(g, ax, pp$t, length(g$widths) - 1, pp$b)

ia2 <- which(g2$layout$name == "ylab-l")
ga2 <- g2$grobs[[ia2]]
ga2$rot <- 90
g <- gtable_add_cols(g, g2$widths[g2$layout[ia2, ]$l], length(g$widths) - 1)
g <- gtable_add_grob(g, ga2, pp$t, length(g$widths) - 1, pp$b)

# draw it

pdf("Voteshare_inc_pm2_total_boot_quad1_firststage_ms2.pdf") # Open a new pdf file
grid.draw(g) # Write the grid.arrange in the file
dev.off() # Close the file


############################
#Figure A18 (c)
############################


fit <- lm(demean_inc_nonpm ~ demean_prev_monsoon2 + I(demean_prev_monsoon2^2), data = Voteshare_inc_nonpm1)
prd <- data.frame(demean_prev_monsoon2 = seq(from = range(Voteshare_inc_nonpm1$demean_prev_monsoon2)[1], to = range(Voteshare_inc_nonpm1$demean_prev_monsoon2)[2], length.out = 100))
err <- predict(fit, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit

quad_Voteshare_inc_nonpm1<-ggplot(prd, aes(x = demean_prev_monsoon2, y = fit)) +
  theme_bw() +
  geom_line() +
  geom_smooth(aes(ymin = lci, ymax = uci), stat = "identity",colour = "black" ) +
  coord_fixed(ratio=0.1) +
  scale_y_continuous(limits = c(-9, 10),breaks = c(-8,-6,-4,-2,0,2,4,6,8,10) ) + 
  #geom_line(aes(y = lo.Voteshare_inc_pm$up.lim)) + geom_line(aes(y = lo.Voteshare_inc_pm$low.lim)) +
  geom_ribbon(data=prd, aes(ymin=lci ,ymax=uci), fill="grey50", alpha="0.5") +
  labs(y = "Vote Share (PM-Incumbent, demeaned)", x="Prev. Monsoon (demeaned)", title="Quadratic Prediction Plot: Vote Share (Non-PM-Incumbent, PM-MP)") + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.title.y=element_text(vjust=-1),
        axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #panel.border = element_blank(),
        #panel.background = element_blank(),
        plot.margin = unit(c(0.3,0.7,0,0), "cm")) 

gg_Voteshare_inc_nonpm1<-ggplot(data=Voteshare_inc_nonpm1, aes(demean_prev_monsoon2)) +  
  geom_histogram(aes(y=..count../sum(..count..)), colour = "grey50", fill = "grey50", breaks=seq(-1, 1, by=0.2), alpha = "0.2") + 
  labs(y = "Share of Observations") +  
  coord_fixed(ratio=1) +
  scale_y_continuous(limits = c(0, 0.4), breaks = c(0,0.1,0.2,0.3,0.4)) + 
  theme(axis.title.y=element_text(vjust=-0.5),
        axis.text.y = element_text(vjust=0),
        #panel.background = element_rect(fill = NA),
        axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.ticks.length=unit(0.2,"cm"),
        plot.margin = unit(c(0.3,0.7,0,0), "cm")) 

grid.draw(gg_Voteshare_inc_nonpm1)

g1 <- ggplot_gtable(ggplot_build(quad_Voteshare_inc_nonpm1))
gpm1 <- ggplot_gtable(ggplot_build(quad_Voteshare_inc_nonpm1))
g2 <- ggplot_gtable(ggplot_build(gg_Voteshare_inc_nonpm1))

pp <- c(subset(g1$layout, name == "panel", se = t:r))
g <- gtable_add_grob(g1, g2$grobs[[which(g2$layout$name == "panel")]], pp$t, 
                     pp$l, pp$b, pp$l)
ia <- which(g2$layout$name == "axis-l")
ga <- g2$grobs[[ia]]
ax <- ga$children[[2]]
ax$widths <- rev(ax$widths)
ax$grobs <- rev(ax$grobs)
ax$grobs[[1]]$x <- ax$grobs[[1]]$x - unit(1, "npc") + unit(0.15, "cm")
g <- gtable_add_cols(g, g2$widths[g2$layout[ia, ]$l], length(g$widths) - 1)
g <- gtable_add_grob(g, ax, pp$t, length(g$widths) - 1, pp$b)

ia2 <- which(g2$layout$name == "ylab-l")
ga2 <- g2$grobs[[ia2]]
ga2$rot <- 90
g <- gtable_add_cols(g, g2$widths[g2$layout[ia2, ]$l], length(g$widths) - 1)
g <- gtable_add_grob(g, ga2, pp$t, length(g$widths) - 1, pp$b)

# draw it

grid.draw(g)

pdf("Voteshare_inc_nonpm1_total_boot_quad1_firststage_ms2.pdf") # Open a new pdf file
grid.draw(g) # Write the grid.arrange in the file
dev.off() # Close the file


############################
#Figure A18 (d)
############################


fit <- lm(demean_inc_nonpm ~ demean_prev_monsoon2 + I(demean_prev_monsoon2^2), data = Voteshare_inc_nonpm0)
prd <- data.frame(demean_prev_monsoon2 = seq(from = range(Voteshare_inc_nonpm0$demean_prev_monsoon2)[1], to = range(Voteshare_inc_nonpm0$demean_prev_monsoon2)[2], length.out = 100))
err <- predict(fit, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit

quad_Voteshare_inc_nonpm2<-ggplot(prd, aes(x = demean_prev_monsoon2, y = fit)) +
  theme_bw() +
  geom_line() +
  geom_smooth(aes(ymin = lci, ymax = uci), stat = "identity",colour = "black" ) +
  coord_fixed(ratio=0.1) +
  scale_y_continuous(limits = c(-9, 10),breaks = c(-8,-6,-4,-2,0,2,4,6,8,10) ) + 
  #geom_line(aes(y = lo.Voteshare_inc_pm$up.lim)) + geom_line(aes(y = lo.Voteshare_inc_pm$low.lim)) +
  geom_ribbon(data=prd, aes(ymin=lci ,ymax=uci), fill="grey50", alpha="0.5") +
  labs(y = "Vote Share (PM-Incumbent, demeaned)", x="Prev. Monsoon (demeaned)", title="Quadratic Prediction Plot: Vote Share (Non-PM-Incumbent, non-PM-MP)") + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.title.y=element_text(vjust=-1),
        axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #panel.border = element_blank(),
        #panel.background = element_blank(),
        plot.margin = unit(c(0.3,0.7,0,0), "cm")) 

gg_Voteshare_inc_nonpm2<-ggplot(data=Voteshare_inc_nonpm0, aes(demean_prev_monsoon2)) +  
  geom_histogram(aes(y=..count../sum(..count..)), colour = "grey50", fill = "grey50", breaks=seq(-1, 1, by=0.2), alpha = "0.2") + 
  labs(y = "Share of Observations") +  
  coord_fixed(ratio=0.1) +
  scale_y_continuous(limits = c(0, 0.4), breaks = c(0,0.1,0.2,0.3,0.4)) + 
  theme(axis.title.y=element_text(vjust=0.5),
        axis.text.y = element_text(vjust=0),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        panel.background = element_rect(fill = NA),
        #axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        #panel.background = element_blank(),
        axis.ticks.length=unit(0.2,"cm"),
        plot.margin = unit(c(0.3,0.7,0,0), "cm")) 



g1 <- ggplot_gtable(ggplot_build(quad_Voteshare_inc_nonpm2))
gpm2 <- ggplot_gtable(ggplot_build(quad_Voteshare_inc_nonpm2))
g2 <- ggplot_gtable(ggplot_build(gg_Voteshare_inc_nonpm2))

pp <- c(subset(g1$layout, name == "panel", se = t:r))
g <- gtable_add_grob(g1, g2$grobs[[which(g2$layout$name == "panel")]], pp$t, 
                     pp$l, pp$b, pp$l)

ia <- which(g2$layout$name == "axis-l")
ga <- g2$grobs[[ia]]
ax <- ga$children[[2]]
ax$widths <- rev(ax$widths)
ax$grobs <- rev(ax$grobs)
ax$grobs[[1]]$x <- ax$grobs[[1]]$x - unit(1, "npc") + unit(0.15, "cm")
g <- gtable_add_cols(g, g2$widths[g2$layout[ia, ]$l], length(g$widths) - 1)
g <- gtable_add_grob(g, ax, pp$t, length(g$widths) - 1, pp$b)

ia2 <- which(g2$layout$name == "ylab-l")
ga2 <- g2$grobs[[ia2]]
ga2$rot <- 90
g <- gtable_add_cols(g, g2$widths[g2$layout[ia2, ]$l], length(g$widths) - 1)
g <- gtable_add_grob(g, ga2, pp$t, length(g$widths) - 1, pp$b)

# draw it

grid.draw(g)

pdf("Voteshare_inc_nonpm2_total_boot_quad1_firststage_ms2.pdf") # Open a new pdf file
grid.draw(g) # Write the grid.arrange in the file
dev.off() # Close the file

